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Abstract - We present a density functional based closure of the pair Smoluchowski equation 
for Brownian particles under shear flow. Given an equilibrium free energy functional as input 
the theory provides first-principles predictions for the flow- distorted pair correlation function and 
associated rheological quantities over a wide range of volume fractions and flow rates. Taking two- 
dimensional hard-disks under shear flow as an illustrative model we calculate the pair correlation 
function, viscosity and normal stress difference under both steady and start-up shear. 



Introduction. — The addition of colloidal particles to 
a Newtonian liquid gives rise to a nonlinear rheological re- 
sponse, characterized by a rate dependent viscosity, finite 
normal stress differences and nontrivial transient dynam- 
ics [l] . Understanding the interplay between particle inter- 
actions and external stress or strain fields remains a major 
theoretical challenge and much effort has been invested in 
the search for tractable closure relations which capture the 
essential physics of systems driven out-of-equilibrium [2]. 
Realistic models for which the competing effects of Brow- 
nian motion, potential and hydrodynamic interactions are 
simultaneously active pose particular difficulties |3j. 

The microstructural distortion induced by an external 
flow field is encoded in the nonequilibrium n-particle cor- 
relation functions. For pairwise interacting Brownian par- 
ticles knowledge of the pair distribution function alone 
is sufficient to calculate the full stress tensor, indicating 
that specific features in this quantity can be correlated 
with nonlinearities in the macroscopic rheology. Such an 
approach was employed in recent experiments by Cheng 
et al. 4 and Koumakis et al. 5 in which confocal mi- 
croscopy was used to analyze the microstructural changes 
associated with the onset of shear thinning, thickening and 
yielding. Brownian and Stokesian dynamics simulations of 
hard spheres have provided further insight [6], but a gen- 
eral theoretical method with firm thermodynamic founda- 
tion is still lacking. Existing theories focus either on states 
close to the glass transition by employing mode-coupling 
approximations [7 10 , or attempt to extend exact low vol- 
ume fraction results lT|[l2 to finite volume fraction using 



liquid state integral equation closures of the pair Smolu- 
chowski equation [13 17 . 

In this paper we present a conceptually simple method 
by which the pair correlations and nonlinear rheology of 
Brownian suspensions can be predicted using dynamical 
density functional theory (DDFT) [l8 19 . The accurate 
treatment of packing effects and robust mean-field descrip- 
tion of phase equilibria provided by modern density func- 
tional approximations are inherited by our theory, thus 
providing a simple and flexible framework within which 
the interplay between particle interactions, diffusion and 
external driving may be studied. An advantage of our 
approach over existing integral equation based theories 
is that we provide a clear link between the macroscopic 
rheology and an underlying equilibrium free energy func- 
tional. The use of an explicit generating function avoids 
the familiar problem of thermodynamic inconsistency and 
no-solution regions of parameter space presented by in- 
tegral equation closures 20 . This feature is of partic- 



ular importance for systems exhibiting equilibrium phase 
transitions, as the proximity of the chosen thermodynamic 
state point to underlying phase boundaries may influence 
the rheological response. The present method is concep- 
tually straightforward and may be applied to calculate the 
viscosity and normal stresses of any model for which there 
exists an explicit approximation of the free energy func- 
tional. It thus becomes possible to exploit the vast array 
of available equilibrium density functional approximations 
(for representative examples see [2l}|24| ) for rheological 
studies. 
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Theory. — We consider a system of N Brownian par- 
ticles interacting via isotropic pairwise interactions, homo- 
geneously dispersed in an incompressible Newtonian fluid 
of given viscosity In the absence of hydrodynamic inter- 
actions the probability distribution of particle positions is 
given by the Smoluchowski equation [3] 
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where Do is the bare diffusion coefficient, f3 = (/c^T) -1 
and the conservative force on particle i is generated from 
the potential energy according to = — V* 5^ ■ u{rij). To 
preserve translational invariance we omit external poten- 
tials and consider a velocity v^(t) = n(t) • r^, where n(t) 
is the traceless velocity gradient tensor. 

Integration of ([!]) over the centre-of-mass of a pair of 
particles and the remaining TV— 2 coordinates leads to the 
pair Smoluchowski equation for the flow distorted pair dis- 
tribution in which only the relative coordinate 1*12 = r 2 — 1*1 
appears 



dg(r 12 ) 
dt 



V-j(r 12 ) = 0, 



(2) 



where V = V2 = — Vi and time arguments have been 
suppressed. The pair flux consists of terms due to affine 
flow, Brownian motion, direct and indirect interactions 

j(ri 2 ) = v(ri 2 )#(ri 2 ) - T k B TVg(r 12 ) + g(r 12 )Vu(r 12 ) 



°-J dr 3 g^\r u r 2 ,r 3 )(\/ 2 u(r 23 ) - Vitz(ri 3 )) 



(3) 



with T = 2j3 Do and p = N/V. The integral terms repre- 
sent the residual influence of the particles which have been 
integrated out and require knowledge of g^ (ri, r 2 , r 3 ). In 
the absence of flow the pair flux vanishes and ([3| reduces 
to the exact second member of the YBG hierarchy [25] . 

The thermodynamic stress can be calculated directly 
from g(r) by integration [26] 



1 f rr 

<r(t) = -k B Tpl + -p 2 J dr-u'(r)g(r,t), 



(4) 



where we have made explicit the fact that all time- 
dependence comes from the pair distribution. It has been 
shown E] for low and intermediate rates that supplement- 
ing the expression Q by a volume fraction dependent high 
frequency viscosity can account for the measured macro- 
scopic viscosity of hard sphere suspensions Q 

In order to close the theory we first note that the triplet 
distribution of the considered homogeneous system may be 
related to an inhomogeneous pair density [27] 



^ (3) (r!,r 2 ,r 3 ) 



(5) 



1 Although the non-hydrodynamic expression Q was employed 
in ^] to calculate the stress, hydrodynamic effects are nevertheless 
present in the experimentally measured g(r) used as input. 
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Fig. 1: Comparison of pair correlation function g(r) between 
simulation (top row) and DDFT (bottom row) at <j) — 0.6. 
The Peclet number employed are Pe — 0.5, 2 and 5, from left 
to right. The qualitative features are very well reproduced, 
although it is evident that the theory slightly underestimates 
the extent of the structural distortion. 



where the subscript indicates that the source of inhomo- 
geneity is the particle located at ri, now viewed formally 
as an external field perturbing the density distribution. 

An equilibrium sum rule may now be used to relate the 
r.h.s of ([5| to the gradient of the one-body direct correla- 
tion function (191. 



^T^(r 12 )V 2 4 1 i )(r 2 ) = Jdr 3 



Pri } (r 2 ,r 3 ) 



V 2 *i(r 23 ), (6) 



where c^(r 2 ) is the direct correlation function at r 2 in 
an external field generated by the particle at ri. The as- 
sumption that (|6| holds also in nonequilibrium constitutes 
an adiabatic approximation [28] , namely that the inhomo- 
geneous pair density relaxes instantaneously to that of an 
equilibrium system with density profile pg(Y\ 2l t). Given 
the Helmholtz free energy 

F[p{v,t)} = k B T J drp(r,t)[ln(A 3 p(r,t)) - 1] 

+ J; x [p(r,t)] + Jdru(r)p(r,t), (7) 

with thermal wavelength A, the direct correlation function 
is generated by differentiation of the excess contribution 



c (D (r) = _^x[p(r)] 



5p(r) 



(8) 



Employing (|5j|8j) and exploiting the symmetry of the direct 
correlation function enables us to approximate the integral 
term in ([3|, yielding a closed equation of motion 



dg(r,t) 
dt 



-V-j(r,t) 



(9) 



j(r,t) = v(r,t)g(r,t)-r g (r,t)V ^ r ' t)] (10) 

dp(r,t) 
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Equations ([7]), (J9| and (10) thus provide a closed theory 
for calculating the pair correlation function which requires 
only that the equilibrium excess free energy functional cor- 
responding to the pair potential u{r) is known. Given 
#(r, £) the stress tensor, and thus the macroscopic rheol- 
ogy, can be obtained from Q. 

In order to test our approach we consider a two dimen- 
sional system of hard disks of radius R subject to shear 
flow v(r, t) = j(t)ye x . This can be considered as a rather 
demanding test case, as the discontinuous nature of the 
interaction potential leads to strong spatial variations in 
g(r) which require both an accurate equilibrium functional 
as input to the dynamical theory and precise numerical al- 
gorithms with sufficient resolution. We employ a recently 
proposed approximation for T ex due to Roth et al. [24] , 
which ensures that the g eq (r) generated by our theory in 
the absence of flow is in excellent agreement with simu- 
lation data. For steady flows the relative importance of 
advection with respect to Brownian motion is quantified 
by the Peclet number Pe = ^R 2 /2Dq, which emerges nat- 
urally from ([9| upon scaling time with 2Do and taking the 
disk radius as the unit of length. Working in two dimen- 
sions enables accurate numerical solution of ([9| using only 
modest computational resources and yields for fluid states 
a phenomenology very similar to that observed in three 
dimensions. 

Numerics. — The numerical treatment of the hard 
disk case poses several difficulties. Firstly, we are primar- 
ily interested in the value of g(r) at contact, i.e. on a circle 
of radius 2R around the center of the fixed test particle 
(see Eq.Q). The requirement to have grid points located 
precisely on discontinuities in g(r) is incompatible with 
uniform cartesian grids of the type usually employed in 
DDFT calculations. Implementation of a cartesian grid 
would require the use of a very fine grid spacing to obtain 
results of sufficient accuracy, resulting in high computa- 
tional demands. Secondly, the no-flux boundary condi- 
tions imposed by the hard-disk test particle require spe- 
cial treatment to ensure that the continuity equation ([9| 
is respected; particle number should not drift significantly 
during the time-evolution as a result of numerical errors. 
Thirdly, the equation is stiff when fine grids are employed, 
so care is required for the time integration. 

To overcome the geometrical incompatibilities of spher- 
ical particles and cartesian grids, we choose a finite- 
element-like discretisation, allowing great flexibility, as the 
grid can be partially refined in regions of interest, such as 
the region of rapid variation in g(r,t) close to contact. 

What distinguishes the present problem from those oc- 
curing in standard computational fluid dynamics is that 
the fundamental measure free energy functional employed 
requires the evaluation of nonlocal spatial convolutions 
which are challenging to evaluate when working on an un- 
structured grid: Fast-Fourier- Transform techniques can- 
not be used. An advantage of the FMT approach is that 
the weight functions are density independent and have lim- 




Fig. 2: Zero-shear viscosity (full line) in units of ksT /2D as 
a function of volume fraction. Circles are simulation data for a 
monodisperse system and crosses are simulation data for a non- 
crystallizing binary mixture (see text). The dot dashed line in- 
dicates an estimate for the random close packing of monodis- 
perse disks (0 rcp = 0.82) [30] . The fit to the monodisperse 
simulation data is given by rj = O.1O90 2 (0.700 - 0)" 917 . In- 
set: The shear rate dependent viscosity as a function of Peclet 
number for volume fractions <j) — 0.1 — 0.5, in steps of 0.1. 



ited range. It is thus possible to precalculate the convolu- 
tion, so that it can be evaluated at each iteration with only 
a matrix- vector product and the solution of a sparse linear 
system, operations that are highly optimized in modern 
numerical linear algebra libraries. We have implemented 
this scheme using the finite element framework deal. II [35] . 

To realise the hard walls both the density and the den- 
sity flux have to be discretized. This enables accurate 
implementation of the boundary condition that the flux 
normal to the hard boundary must vanish, but increases 
the computational effort when compared to a discretiza- 
tion of the density alone. The resulting differential equa- 
tion for p(r, t) is stiff, i.e. its Jacobian has eigenvalues 
with very large eigenvalues stemming from the diffusion 
term 31 . This results in stability issues of the numerical 
solution when using simple timestepping methods, unless 
very small step sizes are chosen. The finer the grid is 
chosen, the more pronounced the stiffness and the smaller 
the steps necessary for stable time integration, dramat- 
ically increasing the computational demand. Often im- 
plicit Runge-Kutta Methods allow a efficient solution of 
stiff differential equations, but in this case the large di- 
mension of the problem and its non-linearity prohibit the 
use of these methods. We thus make use of a class of ex- 
plicit Runge-Kutta Methods that are specifically tailored 
for stiff equations [32] . 

Simulations. — The simulation data were obtained 
using an event-driven algorithm whose detailed descrip- 
tion for the non-sheared case can be found in [33]. Its 
extension to treat external fields is discussed in [34]. We 
simulate a system of N = 1000 hard disks with radius 
R undergoing Brownian motion with the dimensionless 
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Fig. 3: The pair distribution function for Pe = 0.5 (top 
left) and Pe = 5 (top right) at volume fraction <j> — 0.4. We 
also show the angular dependence of fj}~ = sin(0) cos(0)(<7+ — 
g^)/Pe. The viscosity for a given shear rate is given by the 
area under the corresponding curve. The angle is measured 
clockwise, starting from the top. 



diffusion coefficient Dq/(2vqR) = 0.005. Here is gen- 
erated from a Gaussian velocity distribution imposed by 
the thermostat, ensuring that the distribution has the di- 
mensionless variance of moVo/(/c#T) = 2. The Brownian 
time as defined in 33 is set to TBVo/d = 0.01. Equilibra- 
tion runs where performed, ensuring jt ^> 1 (cf. [34]) or 
D t/R 2 > 10 3 (cf. [33 ) respectively. 

In the simulation the pair distribution functions can be 
directly extracted from the particle positions. For the zero 
shear viscosity we use the collision based method from [34] 
which avoids the cumbersome numerical evaluation of the 
limit 7—^0. The viscosity can hence be calculated via 




2k B TV t-> 



(11) 



where r l J stands for the ^/-component of the relative co- 
ordinate between particle i and j and v l J for the x- 
component of the particles' relative velocity, while the sum 
runs over all collisions during the simulation time in the 
simulation box with the volume V. 

Steady shear. — In Fig.l we compare g(r) from 
DDFT with Brownian dynamics data for a volume frac- 
tion (j) = 0.6 and three values of Pe. The chosen vol- 
ume fraction lies well below the crystallization value of 
(j) cr = 0.69 but still represents a dense, strongly interact- 
ing colloidal liquid state. At the lowest shear rate con- 
sidered, Pe = 0.5, the microstructure is only slightly dis- 
torted away from equilibrium and the ring structure of 
nearest, next-nearest etc. neighbour peaks is still clearly 
visible, reflecting the fact that Brownian forces are dom- 
inant over shear forces in this regime. As Pe = 0.5 is 



a close-to-equilibrium state the good level of agreement 
between theoory and simulation owes much to the high 
quality of the excess free energy functional employed [24] , 
which generates the equilibrium structure. The most not- 
icable discrepancy is an overestimation of the higher order 
peak heights in the two extensional quadrants (xy > 0). 
At the intermediate shear rate, Pe = 2, shear forces start 
to dominate the Brownian forces and the microstructure 
becomes qualitatively different from equilibrium. Peaks in 
g(r,t) which lie in the compressional quadrant increase in 
amplitude and are narrowed as particles begin to pile-up 
against the test particle. In contrast, the microstructure 
in the extensional quadrants becomes weaker due to parti- 
cle depletion in the wake behind the test particle. Even at 
this relatively modest value of Pe we observe the forma- 
tion of lane- type structures in the extensional quadrants, 
arising from excluded volume packing constraints. At the 
highest shear rate considered, Pe = 5, the microstruc- 
ture is quite different from that in equilibrium with a very 
pronounced laning structure in the extensional quadrants. 
While the qualitative level of agreement between theory 
and simulation remains good, the packing structure ap- 
pears to be underestimated by the theory, an effect which 
we attribute to the adiabatic approximation used to close 
equation 

In Fig. 2 we show the zero-shear limit of the viscosity 
770 = linLy^o Vxyl - 1 as a function of volume fraction. For 
(j) < 0.45 the calculated 770 is in good quantitative agree- 
ment with Brownian dynamics, but for larger values of (j) 
the simulation data exhibits a rapid increase not captured 
by theory. This discrepancy is due to the onset of slow 
structural relaxation at high volume fractions which be- 
comes partially lost upon making an adiabatic closure of 
the theory. For volume fractions in excess of (j) cr = 0.69 
crystallization becomes an issue. Although our theory is 
in principle sensitive to crystallisation (by virtue of the 
generating functional employed [24 ) the results presented 
here follow the disordered branch of the free energy, which 
in reality becomes metastable above <j) cr . For large values 



of (j) approaching random close packing (0 r , 



0.82 for 



disks 30 ) accurate numerical solution of ([9| becomes dif- 
ficult, but the available data suggests that 770 diverges at 
a volume fraction considerably lower than unity, the value 
at which the input free energy functional becomes singu- 
lar. For comparison we also include in Fig. 2 simulation 
data generated for a binary disc mixture in which half of 
the particles have diameter unity and the other half di- 



ameter 1.4. This mixture is known not to crystallize 34 



and so the close agreement with the monodisperse simula- 
tions suggests that neglect of crystallization in our theory 
is not responsible for the apparent discrepancy. The inset 
to Fig. 2 shows the rate dependent viscosity from DDFT 
for several values of the volume fraction within the range 
< (j> < (j) cr . In each case the Newtonian plateau is fol- 
lowed by a regime of shear thinning characterized by a 
volume fraction dependent thinning exponent. 

In order to investigate the features of g(r) responsible 
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Fig. 4: Stress curves for Start-up shear for <j> — 0.45 and 
different Pe. Stresses are scaled by long time limit. For low 
shear rates the stress increases monotonically towards the long 
time limit. However for Pe larger than a critical value there is 
an overshoot at short times whose relative magnitude increases 
with the shear rate. Inset: Relative magnitude of the overshoot 
for different packing fractions and shear rates. The overshoot 
increases with shear rate and decreases with volume fraction. 



for the shear thinning behaviour we plot in Fig. 3 g(r) at 
a fixed volume fraction, <fi = 0.4, for a Peclet number in 
the linear regime (Pe = 0.5, left panel in the figure) and 
another in the shear thinning regime (Pe = 5, right panel 
in the figure). In determining the rheology of hard particle 
systems only the angular dependence of the contact value 
gt around the particle surface is significant (see equation 
Q). However, closer inspection reveals that it is the quan- 
tity fg~ = sin(0) cos(#)(g+ — g^)/Pe, which contributes to 
the viscosity in Q. In the lower panel of Fig. 3 we show 
/ + as a function of (where the angle is measured in a 
clockwise direction, starting from the vertical). The an- 
gular variation of fg~ clearly shows the reduction in the 
height of the depletion peak at = 7r/4 which, when inte- 
grated, leads to shear thinning. The density aggregation 
in the compressional quadrant at = 37r/4 has only a very 
minor effect on the rheology. 

Startup shear. — In Fig. 4 we show the shear stress 
following the switch on of steady shear, j(t) = 76 (£). This 
time-dependent shear field generates a transient response 
as the system evolves from equilibrium to the steady-state 



36 



For small shear rates the stress monotonically ap- 
proaches its long time value, but for more rapid shearing it 
first increases to a maximum value before slowly decaying 
to the steady state. The stress overshoot is a well known 
effect related to the passage from a regime of primarily 
elastic response to one of dissipative steady-state flow (see 
[5] and 37 for recent investigations of this phenomena). 



Fig. 5: Buildup of normal stress differences (top row) and 
shear modulus G(t) obtained from start-up (bottom row) for 
various shear-rates at cj) — 0.20 (left column) and for various 
packing fractions at Pe = 0.1 (right column). 



Both numerical simulations and experiment have found 
that the height of the overshoot maximum, relative to the 
steady-state shear stress, increases as a function of shear- 



rate but decreases with density [5j[37|- The former effect 
is a rather trivial consequence of appliying a stonger driv- 
ing force, whereas the latter effect arises from the fact 
that suspensions at larger volume fraction have a smaller 
average distance between particle surfaces. Systems foor 
which the particles are densely packed begin plastic flow 
at smaller values of strain and are thus unable to accumu- 
late large elastic stresses prior to the onset of flow. Our 
numerical results, shown in the inset to Fig. 4, are entirely 
consistent with this physical picture. Also consistent with 
the findings of Koumakis et al. is the fact that the criti- 
cal Peclet number for which the overshoot first occurs is 
density dependent (see inset Fig. 4). 

In the upper two panels of Fig. 5 we show the evolution 
of the first normal stress difference Ni = a xx — a yy fol- 
lowing start-up shear flow for various shear rates at fixed 
volume fraction and for various volume fractions at fixed 
shear-rate. The general trends follow closely those of the 
shear stress, but the overshoot is greatly suppressed and is 
not visible on the scale of the figure for the shear rates con- 
sidered. The shear stress following switch-on shear also al- 
lows to obtain the shear modulus simply by differentiation 
with respect to time [2] . It can be seen that the relaxation 
time of the modulus decreases both with increasing shear 
and increasing packing fraction. The bottom right panel 
shows the modulus for <j) = 0.2. We note that at this vol- 
ume fraction the shear modulus, G(t), becomes negative 
at long times for Pe > 1 as a result of the stress overshoot 
in the stress. 

Summary. — In summary, we have developed a test- 
particle approach to calculating the shear-distorted pair 
correlation function of a colloidal suspension. Integration 
of the anisotropy generates the stress tensor, from which 
the nonlinear viscosity and normal stress differences may 
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be calculated without further approximation. The second 
normal stress difference N2 = cr xx — <j yy was not consid- 
ered explicitly here, due to our choice to focus on a two- 
dimensional model system, but is preducted by the theory 
when applied to three-dimensional systems. Furthermore, 
analysis of our fundamental equations, Q, (|9| and (10), 
shows that in three dimensional calculations there will be 
a nonvanishing distortion of g(r) in the vorticity direc- 
tion. Capturing this nontrivial distortion of g(r) has been 
highlighted as an important challenge for theories of the 
nonequilibrium pair correlation functions [16] . Due to our 
rather advanced algorithms for solving the DDFT equa- 
tions we believe that full three dimensional calculations 
are entirely feasible and numerical studies of more realis- 
tic model fluids are the subject of further investigation. 
The approach taken in the present work is close in spirit 



to the micro-rheological studies performed in 38 and 39 



While these works were focused on the density distribution 
about obstacles dragged through a suspension, we have 
focused on using the test particle method to obtain the 
bulk distorted pair-correlations as a gateway to macro- 
scopic rheological studies. A key feature of our work is 
that the main rheological quantities are explicitly con- 
nected to an underlying free energy functional and that 
the method can be generally applied to investigate the 
rheology of any model system for which a free energy ex- 
pression is available. The structural predictions of the 
theory for the difficult test case of hard discs are good 
when compared to the results of simulation (c.f. Fig. 1), 
however, discrepancies are clear in the zero-shear viscosity 
at intermediate and high volume fraction. For hard-discs 
this discrepancy can be attributed to the adiabatic closure 
approximation. By failing to recognize the existence of a 
glass-transition singularity at high density, the theory pre- 
dicts too weak an increase of the viscosity as a function of 
volume fraction. It is interesting that even for monodis- 
perse discs the glass transition has a dominant influence 
on the viscosity for volume fractions below freezing. Given 
that now there exist DDFT-type theories which promise 



to go beyond adiabaticity 40 -42 there would seem to be 



potential to similarly improve upon the theory presented 
here. Work along these lines is ongoing. 
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